Basketball: how likely is it to score?

7.1 μs

xxxxxxxxxx
 
begin
    using Random
    Random.seed!(1234)
end;

10.3 ms

When playing basketball we can ask ourself: how likely is it to score given a position in the court? To answer this question we are going to use data from NBA games from the season 2006 - 2007. We will consider all types of shots.

9.8 μs

xxxxxxxxxx
 
begin
    using CSV
    using DataFrames
    using StatsPlots
    
    season_2007 = CSV.read("./data/seasons/shots_season_2006_2007.csv", DataFrame)
​
    #shots = vcat(season_2007, season_2008, season_2009);
    shots = season_2007
​
    shots[!, :distance] = sqrt.( shots.x .^ 2 +  shots.y .^ 2)
    shots[!, :angle] = atan.( shots.y ./ shots.x )
    filter!(x-> x.distance>1, shots)
end;

37.2 s

DataFrameRow (7 columns)

	result	x	y	period	time	distance	angle
	Int64	Int64	Float64	Int64	Time…	Float64	Float64
1	0	-14	16.75	1	11:42:00	21.8303	-0.874592

xxxxxxxxxx
 
first(shots)

55.1 μs

But how do we interpret the data?

Below we show a sketch of a basketball court, its dimensions and how to interpret the data in the table.

6.7 μs

xxxxxxxxxx
 
begin
    using Images
    img = load("./images/basket_court.png")
end

7.9 s

So, the x and y axis have their origin at the hoop, and we compute the distance from this point to where the shot was made. Also, we compute the angle with respect to the x axis, showed as θ in the sketch. In the data we have the period, which can take values from 1 to 4, meaning the period in which the shot was made.

1.5 ms

We now plot where the shots where made:

5.2 μs

xxxxxxxxxx
 
begin
    histogram2d(shots.y[1:10000], shots.x[1:10000], bins=(50,30))
    ylabel!("y axis")
    xlabel!("x axis")
end

9.7 s

We see that the shots are very uniformly distributed near the hoop, except for distances very near to the hoop, to see this better, we plot the histograms for each axis, x and y.

6.8 μs

But we are interested in the shots that were scored, so we filter now the shots made and plot the histogram of each axis.

2.6 μs

xxxxxxxxxx
 
shots_made = filter(x->x.result==1, shots);

312 ms

xxxxxxxxxx
 
begin
    histogram(shots_made.y[1:10000], legend=false, nbins=40)
    xlabel!("x axis")
    ylabel!("Counts")
end

758 ms

xxxxxxxxxx
 
begin
    histogram(shots_made.x[1:10000], legend=false, nbins=45)
    xlabel!("y axis")
    ylabel!("Counts")
end

21.9 ms

If we plot a 3d plot of the count we obtain the plot wireplot shown below.

2.9 μs

xxxxxxxxxx
 
begin
    using StatsBase
    h = fit(Histogram, (shots_made.y, shots_made.x), nbins=40)
    #plot(h) # same as histogram2d
    wireframe(midpoints(h.edges[2]), midpoints(h.edges[1]), h.weights, zlabel="counts", xlabel="y", ylabel="x", camera=(40,40))
title!("Histogram of shots scored")
end

1.1 s

We see that more shot are made as we get near the hoop, as expected.

It is important to notice that we are not showing the probability of scoring, we are just showing the distribution of shot scored, not how likely is it to score.

3.7 μs

Modeling the probability of scoring

2.8 μs

The first model we are going to propose is a Bernoulli model.

Why a Bernoulli Distribution?

A Bernoulli Distribution results from an experiment in which we have 2 possible outcomes, one that we usually called a success and another called a fail. In our case our success is scoring the shot and the other possible event is failing it.

The only parameter needed in a bernoulli distribution is the probability p of having a success. We are going to model this parameter as a logistic function:

13.9 μs

xxxxxxxxxx
 
begin
    using Turing
    using StatsFuns: logistic
end

23.9 s

xxxxxxxxxx
 
begin
plot(logistic, legend=false)
ylabel!("Probability")
xlabel!("x")
title!("Logistic function (x)")
end

1.3 s

Why a logistic function?

We are going to model the probability of shoot as a function of some variables, for example the distance to the hoop, and we want that our probability of scoring increases as we are getting closer to it. Also out probability needs to be between 0 an 1, so a nice function to map our values is the logistic function.

5.9 μs

So, the model we are going to propose is:

$p \sim l o g i s t i c (a + b * d i s t a n c e [i] + c * a n g l e [i])$

$o u t c o m e [i] \sim B e r n o u l l i (p)$

x
md" So, the model we are going to propose is: 
​
$p\sim logistic(a+ b*distance[i] + c*angle[i])$
$outcome[i]\sim Bernoulli(p)$
"

3.2 μs

But what values and prior distributions are we going to propose to the parameters a, b and c?

Let's see:

9.1 μs

Prior Predictive Checks: Part I

3.1 μs

Suppose we say that our prior distributions for a, b and c are going to be 3 gaussian distributions with mean 0 and variance 1. Lets sample and see what are the possible posterior distributions for our probability of scoring p:

$a \sim N (0, 1)$

$b \sim N (0, 1)$

$c \sim N (0, 1)$

11.1 μs

168 ms

xxxxxxxxxx
 
begin
    plot(possible_distances, predicted_p[1], legend = false, color="blue")
    for i in 2:n_samples
        plot!(possible_distances, predicted_p[i], color=:blue)
    end
    xlabel!("Normalized distance")
    ylabel!("Predicted probability")
    title!("Prior predictive values for p")
    current()
end

258 ms

We see that some of the predicted behaviours for p don't make sense. For example, if b takes positive values, we are saying that as we increase our distance from the hoop, the probability of scoring also increase. So we propose instead the parameter b to be the negative values of a LogNormal distribution. The predicted values for p are shown below.

7.5 μs

So our model now have as priors distributions:

$a \sim N o r m a l (0, 1)$

$b \sim L o g N o r m a l (1, 0.25)$

$c \sim N o r m a l (0, 1)$

3.5 μs

and sampling values from those prior distributions, we obtain the plot shown below for the predicted values of p.

7.6 μs

37.4 ms

70.7 ms

Now that we have the expected behaviour for p, we define our model and calculate the posterior distributions with our data points.

7.5 μs

Defining our model and computing posteriors

5.5 μs

Now we define our model to sample from it:

4.2 μs

logistic_regression (generic function with 1 method)

xxxxxxxxxx
 
@model logistic_regression(distances, angles, result,n) = begin
  N = length(distances)
  # Set priors.
    a ~ Normal(0,1) 
    b ~ LogNormal(1,0.25)
    c ~ Normal(0,1)
    for i in 1:n
        p = logistic( a -  b*distances[i] + c*angles[i])
        result[i] ~ Bernoulli(p)
    end
end

94.3 μs

xxxxxxxxxx
 
n=1000;

42.0 ns

The output of the sampling tell us also some information about sampled values for our parameters, like the mean, the standard deviation and some other computations.

5.4 μs

xxxxxxxxxx
 
# Sample using HMC.
chain = mapreduce(c -> sample(logistic_regression(shots.distance[1:n] ./ maximum(shots.distance[1:n] ), shots.angle[1:n], shots.result[1:n], n), NUTS(), 1500),
    chainscat,
    1:3
);
​

74.4 s

Traceplot

In the plot below we show a traceplot of the sampling.

What is a traceplot?

When we run a model and calculate the posterior, we obtain sampled values from the posterior distributions. We can tell our sampler how many sampled values we want. A traceplot is just showing them in sequential order. We also can plot the distribution of those values, and this is what is showed next to each traceplot.

11.4 μs

xxxxxxxxxx
 
plot(chain, dpi=60)

5.0 s

71.6 ms

Now plotting the probability of scoring using the posterior distributions of a, b and c for an angle of 45°, we obtain:

11.5 μs

516 ms

588 ms

The plot shows that the probability of scoring is higher as our distance to the hoop decrease, which makes sense, since the difficulty of scoring increase.

We plot now how the probability varies with the angle for a given distance. Here we plot for a mid distance, corresponding to 0.5 in a normalized distance.

6.5 μs

566 ms

843 ms

We see that the model predict an almost constant probability for the angle.

3.7 μs

New model and prior predictive checks: Part II

2.8 μs

Now we propose another model with the form:

$p \sim l o g i s t i c (a + b^{d i s t a n c e [i]} + c * a n g l e [i])$

*But for what values of b the model makes sense?

We show below the plot for 4 function with 4 possible values of b, having in mind that the values of x, the normalized distance, goes from 0 to 1.

x
 
md"Now we propose another model with the form: 
​
$p\sim logistic(a+ b^{distance[i]} + c*angle[i])$
​
*But for what values of b the model makes sense?
​
We show below the plot for 4 function with 4 possible values of *b*, having in mind that the values of *x*, the normalized distance, goes from 0 to 1. 
"

6.4 μs

xxxxxxxxxx
 
begin
    f1(x) = 0.3^x
    f2(x) = 1.5^x
    f3(x) = -0.3^x
    f4(x) = -1.5^x
end;

45.6 μs

485 ms

Analysing the possible values for b, the one that makes sense is the value proposed in f1, since we want an increasing influence of the distance in the values of p as the distance decreases, since the logistic function has higher values for higher values of x.

5.3 μs

So now that we know the values the our parameter b can take, we propose for it a beta distribution with parameters α=2 and β=5, showed in the plot below.

4.7 μs

518 ms

Defining the new model and computing posteriors

2.6 μs

We define then our model and calculate the posterior as before.

2.4 μs

logistic_regression_exp (generic function with 1 method)

xxxxxxxxxx
 
@model logistic_regression_exp(distances, angles, result, n) = begin
  N = length(distances)
  # Set priors.
    a ~ Normal(0,1) 
    b ~ Beta(2,5)
    c ~ Normal(0,1)
    for i in 1:n
        p = logistic( a +  b .^ distances[i] + c*angles[i])
        result[i] ~ Bernoulli(p)
    end
end

110 μs

xxxxxxxxxx
 
# Sample using HMC.
chain_exp = mapreduce(c -> sample(logistic_regression_exp(shots.distance[1:n] ./ maximum(shots.distance[1:n] ), shots.angle[1:n], shots.result[1:n], n), HMC(0.05, 10), 1500),
    chainscat,
    1:3
);
​

64.2 s

Plotting the traceplot we see again that the variable angle has little importance since the parameter c, that can be related to the importance of the angle variable for the probability of scoring, is centered at 0.

4.9 μs

57.1 ms

572 ms

Employing the posteriors distributions computed, we plot the probability of scoring as function of the normalized distance and obtain the plot shown below.

2.5 μs

479 ms

Given that we have 2 variables, we can plot the mean probability of scoring as function of the two and obtain a surface plot. We show this below.

2.7 μs

449 ms

im3

252 ms

The plot show the behaviour expected, an increasing probability of scoring as we get near the hoop. We also see that there is almost no variation of the probability with the angle.

2.5 μs

Does the Period affect the probability of scoring?

2.8 μs

Now we will try to answer this question. We propose then a model, and calculate the posterior for its parameters with data of one of each of the four possible periods. We define the same model for all four periods. Also, we don't take into account now the angle variable, since we have seen before that this variable is of little importance.

We filter then our data by its period and proceed to estimate our posterior distributions.

4.0 μs

xxxxxxxxxx
 
shots_period1= filter(x->x.period==1, shots);

59.7 ms

logistic_regression_period (generic function with 1 method)

xxxxxxxxxx
 
@model logistic_regression_period(distances, result,n) = begin
  N = length(distances)
  # Set priors.
    a ~ Normal(0,1) 
    b ~ Beta(2,5)
    for i in 1:n
        p = logistic( a +  b .^ distances[i])
        result[i] ~ Bernoulli(p)
            end
end

129 μs

42.0 ns

xxxxxxxxxx
 
# Sample using HMC.
chain_period1 = mapreduce(c -> sample(logistic_regression_period(shots_period1.distance[1:n_] ./ maximum(shots_period1.distance[1:n_] ),shots_period1.result[1:n_], n_), HMC(0.05, 10), 1500),
    chainscat,
    1:3
);

31.2 s

xxxxxxxxxx
 
shots_period2= filter(x->x.period==2, shots);

60.4 ms

xxxxxxxxxx
 
# Sample using HMC.
chain_period2 = mapreduce(c -> sample(logistic_regression_period(shots_period2.distance[1:n_] ./ maximum(shots_period2.distance[1:n_] ), shots_period2.result[1:n_], n_), HMC(0.05, 10), 1500),
    chainscat,
    1:3
);

29.7 s

xxxxxxxxxx
 
shots_period3= filter(x->x.period==3, shots);

65.0 ms

xxxxxxxxxx
 
# Sample using HMC.
chain_period3 = mapreduce(c -> sample(logistic_regression_period(shots_period3.distance[1:n_] ./ maximum(shots_period3.distance[1:n_] ), shots_period3.result[1:n_], n_), HMC(0.05, 10), 1500),
    chainscat,
    1:3
);

29.7 s

xxxxxxxxxx
 
shots_period4 = filter(x->x.period==4, shots);

76.4 ms

xxxxxxxxxx
 
# Sample using HMC.
chain_period4 = mapreduce(c -> sample(logistic_regression_period(shots_period4.distance[1:n_] ./ maximum(shots_period4.distance[1:n_]), shots_period4.result[1:n_], n_), HMC(0.05, 10), 1500),
    chainscat,
    1:3
);

28.8 s

139 ms

We plot now for each period the probability of scoring for each period, each mean and one standard deviation from it.

3.0 μs

162 ms

Finally, we see that for the periods 1 and 4, the first and the last periods, the probabity of scoring is slightly higher than the other two periods, meaning that players are somewhat better scoring in those periods.

4.3 μs

Summary

In this chapter, we used the NBA shooting data of the season 2006-2007 to analyze how the scoring probability is affected by some variables, such as the distance from the hoop and the angle of shooting.

First, we inspected the data plotting a heatplot of all the shots made and histograms of the scored ones. As our goal was to study the scoring probability, a Bernoulli trial situation, we decided to use a Bernoulli model. Since the only parameter needed in a Bernoulli distribution is the probability p of having a success, we modeled p as a logistic function: $p \sim l o g i s t i c (a + b * d i s t a n c e [i] + c * a n g l e [i])$

We set the prior probability of the parameters a and c to a normal distribution and * b* to a log-normal one. So we constructed our logistic regression model and sampled it using the Markov Monte Carlo algorithm. To get a better understanding of the sampling process we made a traceplot that shows the sampled values in a sequential order.

Later, we decided to try with a more complex logistic regression model, similar to the first one but this time modifying the distance parameter: $p \sim l o g i s t i c (a + b^{d i s t a n c e [i]} + c * a n g l e [i])$

We set the prior distribution of b to a beta distribution and constructed the second logistic regression model, sampled it and plotted the results.

Finally, we saw if the scoring percentage is the same in all the periods.

16.4 μs